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We discuss the procedure for the exact solution of the Riemann problem in special relativistic 
magnetohydrodynamics (MHD). We consider both initial states leading to a set of only three 
waves analogous to the ones in relativistic hydrodynamics, as well as generic initial states leading 
to the full set of seven MHD waves. Because of its generality, the solution presented here could 
serve as an important test for those numerical codes solving the MHD equations in relativistic 
regimesf. 



1. Introduction 

As first formulated by Riemann more than a hundred years ago, the solution of the one- 
dimensional Riemann problem in hydrodynamics consists of determining the temporal evolution 
of a fluid which, at some initial time, has two adjacent uniform states characterized by different 
values of uniform velocity, pressure and density. These initial conditions completely determine 
the way in which the discontinuity will decay after removal of the barrier separating the initial 
"left" and "right" states. 

The Riemann problem has ceased to be merely academic and has gained enormous importance 
when it was realized that its numerical solution can serve as the building block of hydrodynamical 
codes based on Godunov-type finite difference methods JGodunov 1959> . In such methods, the 
computational domain is discretized and each interface between two adjacent grid-zones is used 
to construct the initial left and right states of a "local" Riemann problem. The evolution of the 
hydrodynamical equations is then obtained through the solution across the computational grid of 
the sequence of local Riemann problems set up at the interfaces between successive grid-zones 
(see Godunov 1959. but also .Marti & Miiller 2003. and Font 2003 for the use of the Riemann 
problem in relativistic regimes). 

In general, the Riemann problem requires the solution of a nonlinear algebraic system of equa- 
tions written as a function of a single unknown quantity {e.g. the total pressure at the contact 
discontinuity in purely hydrodynamical problems). With the exception of few trivial initial con- 
figurations, the solution of the Riemann problem cannot be obtained analytically but requires a 
numerical approach. The solution found in this way is referred to as the "exact" solution of the 
Riemann problem, to distinguish it from the "approximate " solution of the Riemann problem, 
which is instead obtained when the system of equations is reduced to a locally-linear form (an 
exhaustive discussion of approximate Riemann solvers can be found in lToro 19 99 1. It is therefore 
useful to stress that although named "exact", the solution of the Riemann problem is necessarily 
obtained with a small but nonzero truncation error. 

The exact solution of the Riemann problem in relativistic hydrodynamics has been obtained 

f The numerical code computing the exact solution is available from the authors upon request. 
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only rather recently and was proposed by Marti & M iiller (1994) for flows that are purely along 
the direction normal to the initial discontinuity. This work has then been extended to the case 
in which tangential velocities are present ( Po ns et al. 2000> and improved in efficiency by ex- 
ploiting the relativistic invariant relative velocity between the two states to predict the wave 
pattern produced ( Rezzo lla & Zanotti 2001 and RezzoUa et al. 2003 1. The relevance of these cal- 
culations has not been restricted to fundamental issues of relativistic hydrodynamics. Quite the 
opposite, these solutions have been of great importance for the testing of complex multidimen- 
sional codes implementing High Resolution Shock Capturing (HRSC) methods, and that are 
based on the approximate or exact solution of Riemann problems at the interfaces between the 
numerical cells ( LeVeque 1992V These codes have then been used in various simulations in ei- 
ther fixed (,Aloy et aLJ^99^Fbnt&^^ai£ne20 et al. 2003t or dynamical spacetimes 

dPuez et al. 2004l|Shibata & Sekiguchi 2005|IBaiotti et al. 2005t . 

This intense and recent development of numerical codes for the solution of the relativistic hy- 
drodynamic equations has been accompanied by an equally intense development of codes solving 
the equations of magnetohydrodynamics (MHD) in relativistic regimes. The reason behind this 
activity is the widespread expectation that strong magnetic fields are crucial in the study and ex- 
planation of several puzzling astrophysical phenomena such as relativistic jets or 7-ray bursts. As 
a result, and in the hope of clarifying issues in relativistic astrophysics which cannot be described 
satisfactorily through analytic techniques, several groups have recently constructed numerical 
codes solving the equations of relativistic MHD on either fixed spacetimes (see, for &x?Lm.- 
ple. IDel Zanna et al. 2003l and lKomissarov 1999l for a flat background and lGammie et al. 20031 
IDe Villiers et al. 2003 Komi ssarov 2004llMizuno et al. 20 04 Fragile 2005 Anton et al. 2005l for 
a black-hole background) or in fully dynamical spacetimes (Duez et al. 2005 1. 

Just like their hydrodynamical counterparts, some of these codes are based on the solution 
of a local Riemann problem suitably formulated for a magnetized fluid, and all are meant to 
be used for ultrarelativistic flows. However, unlike their hydrodynamical counterparts, these 
codes cannot benefit from the comparison with the exact solution of the Riemann problem in 
relativistic MHD. The literature on the Riemann problem in MHD is, in fact, much more lim- 
ited and a general exact solution was found rather recently and for a Newtonian fluid only 
l |Ryu & Jones 1995] IFalle eF al. 1998 1. The background knowledge in this area is even more 
scarce for a relativistic fluid and while no general exact solution has been proposed yet, recent 
work has been made to derive an exact solution in the particular case in which the magnetic 
field of the initial states is tangential to the discontinuity and orthogonal to the fluid veloc- 
ity ( Rom ero et al. 20 05 1. Besides having a larger set of equations when compared to the cor- 
responding problem in relativistic hydrodynamics, a considerable addition to the complexity 
of the Riemann problem in relativistic MHD is represented by the fact that the mathematical 
structure of the problem itself is modified and the system of equations is no longer strictly 
hyperbolic (JLichnerowicz 1967t t. The possibility of having coincident eigenvalues poses the 
question of the uniqueness of the solutions and this represents then a problem within the prob- 
lem. As we will comment also later on, a lively debate on these issues is presently ongoing 
and progress is starting to be made, although first results are known in Newtonian MHD only 
(see Torrilhon 2003b 1. Because the focus of this work is the exact solution of the Riemann prob- 
lem in relativistic MHD as an aid to the development of numerical codes, hereafter we will adopt 
the working assumption that the Riemann problems considered here have a solution and that this 
solution is unique. Clearly, this hypothesis avoids the issue rather than solving it, but it allows for 



f We recall that a systems of m quasi-linear partial differential equations is said to be hyperbolic if 
the matrix of coefficients has m real eigenvalues; furthermore, the system is said to be totally or strictly 
hyperbolic if the eigenvalues are real and also all distinct. 
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a marked progress at least in those cases in which compound waves are not found in numerically 
approximate solutions. 

A direct and important consequence of the scarcity of works in this area of fundamental rela- 
tivistic MHD is that the validation of modern complex MHD codes for the most elementary and 
yet demanding tests has not been made in a quantitative manner for generic initial conditions. 
Rather, it has taken place through the qualitative comparison with the large set of test-problems 
in relativistic MHD meticulously collected over the years (see, for instance, 'Komiss arov 19991 
and Balsara 2001 ) . It should be recognized, however, that for non-generic initial states it is suf- 
ficient to have exact solutions for MHD shocks and rarefactions as this covers all types of basic 
hyperbolic waves of the system and indeed exact solutions of this type were used by Komissarov 
(1999) for quantitative testing. 

On the other hand, the purpose of this paper is to present the procedure for the exact solution of 
the Riemann problem in relativistic MHD with generic initial conditions. Our approach considers 
both initial states with a zero component of the magnetic field along the flow and leading to a set 
of only three waves analogous to the ones in relativistic hydrodynamics, as well as generic initial 
states leading to the full set of seven MHD waves. The approach discussed for the numerical 
solution is based on a "hybrid" approach which adopts different sets of equations according to 
the values of the normal magnetic field and that has turned out to be crucial for a successful 
solution. 

The paper is organized as follows: Section|2]contains the basic equations of relativistic MHD, 
while Section|3ldescribes the strategy used to solve the Riemann problem numerically and which 
combines the methods discussed in Sectionl^and Section|5l Section|6lfocuses on the details of 
the numerical implementation and discusses the solution of a number of tests that have become 
standard references. Finally, the conclusions are collected in Sectional 

We use a spacelike signature (—,+,+,+) and a system of units in which c = 1. Greek indices 
are taken to run from to 3, Latin indices from 1 to 3 and we adopt the standard convention for 
the summation over repeated indices. Finally we indicate 3-vectors with an arrow and use bold 
letters to denote 4- vectors and tensors. 



2. Equations of Relativistic MHD 

Consider an ideal but magnetized relativistic fluid with an energy-momentum tensor given by 

T'^"^ - (p + pe + + 2p,„) u^^u'' + (pg + p„0 ry'^'^ - 6^6'' , (2.1) 

where p is the rest mass density, e the specific internal energy, pg the gas pressure, pm the mag- 
netic pressure, = W(l, w^, w^, w^) the four-velocity, W = — v'^Vi = l/\/l — t;^ the 
Lorentz factor and the 4-vector b has components 

b''=i^Wiv-B),^+Wiv-B)v'^. (2.2) 

Here B is the magnetic field 3-vector and 

fe2 = b'b, ^i^ + i^-Bf = 2Pm . (2.3) 

The general relativistic equations of MHD are then simply obtained after requiring the conser- 
vation of baryon number 

V^(pu^) = 0, (2.4) 
where V represents a covariant derivative, the conservation of energy and momentum 

V^T'"' = , (2.5) 



4 



B. Giacomazzo & L. Rezzolla 



together with the relevant pair of Maxwell equations. If the fluid is assumed to have an infinite 
electrical conductivity (i.e. ideal MHD limit), the Maxwell equations reduce to 9[QFg^] = 0, 
where F is the Faraday tensor and the square brackets refer to antisymmetrised indices. Using 
the definition i2.2\ . the Maxwell equations can be simply written as 



(2.6) 



The system of equations is completed with an equation of state (EOS) relating the 

pressure to the rest-mass density and/or to the energy density. Although hereafter we will use an 
ideal-gas EOS: pg = pe{T — 1), where F is the polytropic index, the procedure described for the 
solution of the Riemann problem is valid for a generic EOS. 

We next assume that the system has a planar-symmetry, i.e. that in a Cartesian coordinate 
system {t, x, y, z) all the variables depend only on t and x, and that the spacetime is flat so that 
covariant derivatives in equations ( I2.4t -( lz!6l l can be replaced by partial derivatives and = Ai 
for any 3-vector A. In this case, the complete set of MHD equations can be written as a set of 
first-order partial differential equations in a flux-conservative form 

1^7 + TT = , (2.7) 
at ox 

where U and F are respectively the vectors of conserved quantities and fluxes, defined as 



U = 



D 

T - 606" 

sv _ b°by 

S"- - 6%^ 



\ 



I 



( Dv"" 

svy^ - h^by 

gz^x _ j^x^z 

BVy^ _ B^vy 

QZ^X _ QXyZ 



\ 



(2.8) 



and where the following definitions have been used 
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D 
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pW , 



P 



D 
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': ph , 



Pn 



h = 



= 1 



Pg 
P 



P 



(2.9) 
(2.10) 
(2.11) 
(2.12) 
(2.13) 

(2.14) 



where h is the total specific enthalpy and hg the one of the gas only. 

Note that the divergence-free condition for the magnetic field and the Maxwell equation for 
the evolution of the x-component of the magnetic field imply that dtB^ = = dxB^, i.e. B^ is 
uniform in space, constant in time and thus always maintains its initial values. 



3. Strategy of Solution 

The general Riemann problem in relativistic MHD consists of a set of seven nonlinear waves: 
tv/o fast-waves (FW), two slow-waves (SW), Wo Alfven-waves (AW), and a contact discontinuity 
(CD) at which only the density may be discontinuous. The fast and slow nonlinear waves can be 
either shocks or rarefaction waves, depending on the change in the pressure and in the norm of 
the magnetic field across the wave. 
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Figure 1. Spacetime structure of the MHD Riemann problem in the case in which the magnetic field 
has tangential components only, i.e. = 0. The "Riemann-fan" in this case is composed of only two 
fast-waves (FW) and of a central tangential discontinuity (TD), thus resembling structure of the Riemann 
problem in pure hydrodynamics. Indicated with R1-R4 are the 4 different regions into which the Riemann 
problem can be decomposed, each representing a different state. 



Building on the experience with relativistic hydrodynamics, our general strategy in the search 
for the solution consists of expressing all of the variables after each wave as functions of the 
values of the same variables ahead of the wave and of an unknown variable behind the wave. 
When considering the Riemann problem in relativistic hydrodynamics, in fact, the solution is 
found after expressing all of the quantities behind the wave as functions of the value of the 
pressure at the contact discontinuity. In this way, the problem is reduced to the search for the 
value of the pressure that satisfies the jump conditions at the contact discontinuity. 

When considering the Riemann problem in relativistic MHD, on the other hand, two differ- 
ent cases need to be distinguished. Assuming the initial discontinuity to have normal along the 
a:-axis, the initial magnetic field in this direction can either be zero {i.e. = 0) or not (i.e. 

^ 0). In the first case, the structure of the solution is very similar to the hydrodynamical one, 
with only two fast waves and a tangential discontinuity along which only the total pressure and 
the X component of the velocity are continuous. The spacetime structure of the Riemann problem 
in this case is sketched in Figure [2 where the "Riemann-fan" is shown to be composed of only 
two fast-waves (FW) and of a central tangential discontinuity (TD). Because of this analogy, the 
numerical solution of the Riemann problem when B''^ = follows the same procedure imple- 
mented in relativistic hydrodynamics. We refer to this as the "total-pressure approach" or simply, 
the "p-method" . 

A detailed investigation of the exact solution of the Riemann problem with tangential magnetic 
fields and when the additional condition v ■ B = is imposed, has been recently proposed 
by Romero et al. (2005). Among the many points discussed, this work has shown that when 
B^ — = V ■ B the Riemann problem in relativistic MHD can be assimilated to the one in 
relativistic hydrodynamics and that all of the corrections introduced by the magnetic field can be 
incorporated in the definition of a new, effective EOS. 

In the second case, on the other hand, the Riemann problem is considerably more complex 
and all of the seven waves are allowed to form when the initial discontinuity is removed. The 
spacetime structure of the Riemann problem in this case is sketched in Figure 121 where the 



6 B. Giacomazzo & L. Rezzolla 




Figure 2. Spacetime structure of the MHD Riemann problem in the general case in which the magnetic 
field has also a normal component, i.e. 7^ 0. The "Riemann-fan" is here composed of two fast-waves 
(FW), of two Alfven waves (AW), of two slow-waves (SW) and of a central contact discontinuity (CD). 
Indicated with R1-R8 are the 8 different regions into which the Riemann problem can be decomposed, 
each representing a different state. Indicated are also the different methods used to compute the solutions in 
the different regions (i.e. Bt-method in regions R4 and R5 and p-method in regions R2-R3 and R6-R7). 

"Riemann-fan" is shown to be composed of two fast-waves (FW), of two Alfven waves (AW), of 
two slow-waves (SW) and of a central contact discontinuity (CD). 

It is important to bear in mind that across the Alfven discontinuities only the total pressure, the 
gas pressure and the density are continuous, while there could be jumps in the other quantities. As 
a result, if the total pressure is used as unknown, there would be three different values for the total 
pressure (two between the fast and the slow-waves and one between the two slow-waves) but five 
conditions to be satisfied at the contact discontinuity: the continuity of the three components of 
the velocity and the continuity of the tangential components of the magnetic field. The resulting 
system of five equations in three unknowns is over-constrained and there is no guarantee that a 
global convergent solution is found at the contact discontinuity. Indeed, experience has shown 
that small numerical imprecisions at the level of round-off errors are in general sufficient to 
prevent the simultaneous solution of the five constraints. 

To circumvent this difficulty and inspired by the procedure followed in the exact solution of 
the corresponding Riemann solver in nonrelativistic MHD ( Ryu & Jones 199^ , when i?^ ^ 
we have implemented a "hybrid" approach in which the total pressure is used as the unknown 
variable between the fast and the slow waves {i.e. in regions R2-R3, and R6-R7 of Figure |2}, 
while the tangential components of the magnetic fields (i?^ and B^') is used between the slow 
waves {i.e. in regions R4-R5 of Figure|2jl. In this way, the continuity of the tangential components 
of the magnetic field By and is automatically guaranteed through the contact discontinuity 
and only the continuity of the total pressure and of the three components of the velocity needs 
to be satisfied. The resulting system consists of four equations in four unknowns and, being 
closed, it can be solved numerically through root-finding techniques for nonlinear system of 
equations {e.g. using a Newton-Raphson method). We refer to this as the "tangential magnetic 
field approach" or simply, the "Bt-method" . 

As mentioned in the Introduction, hereafter we will assume that the Riemann problem has a 
solution and that this is unique. As a result, we will not discuss in any detail compound waves 
which seem to develop in the numerical solution of some special initial states (one of these is 
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shown in Section 16.2. U and whose admissibility as solution of the Riemann problem is still 
debated. 



4. Total-Pressure Approach: "p-method" 

In the following Sections we describe in detail the approach in which we calculate all of the 
variables in the Riemann fan using as unknown the total pressure, i.e. the p-method. Different set 
of equations will be derived according to whether the solution is across a shock or a rarefaction 
wave. 



4. 1 . Solution across a shock front 

Consider E to be a hypersurface in flat spacetime across which p, u and T are discontinuous. Let 
also n be the unit 4-vector normal to E so that the Rankine-Hugoniot conditions for relativistic 
MHD can be expressed as 

[[pu"]]n„ = 0, (4.1) 
[[T"'3]]„„ = 0, (4.2) 
[[6"w'3-w"&'']]n„ = 0, (4.3) 

where we use the double-bracket notation to express the jump of a quantity F across the hyper- 
surface S, i.e. 

[[F]] =Fa-Fi,. 

where Fa and Ff, are respectively the values ahead (a) and behind (b) the shock. 

In particular, if S is the 4-dimensional hypersurface describing the evolution of a shock wave 
normal to the a;-axis, the unitary condition on n can be used to derive the components 

n'^ = WsiVs, 1,0,0) , (4.4) 

where Vs is the coordinate velocity of the shock, Ws = (1 — Vg)~^^^ its Lorentz factor, and we 
can rewrite equations (I4.1> -( l43t explicitly as 



[[J]] ^ [[pW{Vs - 


V^Ws]] 


= 0, 


(4.5) 


[[b%°~T]]Vs+[[S''-b%''- 


-Dv^]] 


= 0, 


(4.6) 


[[b°b^ - S'']]Vs + [[S'^v'' + p 




= 0, 


(4.7) 


[[b%y-sy]] Vs + [[syv^ 


- b'^by]] 


= 0, 


(4.8) 






= 0, 


(4.9) 




[[B^]] 


= 0, 


(4.10) 


[[By]] Vs + [[B-yy - 


-v'^By]] 


= 0, 


(4.11) 


[[B^]] Vs + [[B-v' - 




= 0, 


(4.12) 



where J is the (rest) mass flux across the shock. 

After a number of tedious but otherwise straightforward algebraic manipulations, equations 
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I4.5>-(|4TT2> can be recast as 
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w^ 



J 

w. 
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- 




By ■ 
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W^rj^v^ 
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' B^ ' 






D 








D 






w^_ 





D 

B- - ^ 

B^ [[..1] - ^ 



1 



0, (4.13) 

0, (4.14) 

0, (4.15) 

0, (4.16) 

0, (4.17) 



51/ 

U 

[[B-]]=0, (4.18) 



J 




By 






D 


J 




'B^ 


Ws 




D 



B^ \\v'] 



0, (4.19) 
0, (4.20) 



where we have defined rj = v ■ B and exploited the property 



mvs - 



J 

w. 



valid for any scalar quantity F. 

The next step to take is to express all of the variables as functions of J and pi, only. We start 
by using equation (I4.13> to obtain 



SO that equation (I4.14> yields 



n 



W 

— [B"" iVa - Vb) - PaVl + PbVl] 



(4.21) 



(4.22) 



which depends on vf, pb but also on B^, B^, vf, v^. To remove the dependence from these latter 
quantities we employ equations ( I4.19> and (I4.20> to obtain Bj^ and B^ as functions of v^, vf 
and pb, i.e. 



Bl 



Db 



Bl 


w, 

h -^B^vl - 


Ws 


Da 


J 


J 


Bl 


W, 


Ws 


Da 




J 



B^vi 



(4.23) 
(4.24) 



We can now solve equation (I4.15> and finally obtain as a function of v^,v^, pb and J 

Da {BlWs + Wl [Ws{pb - Pa) - BlWsjl - yyaVl - Vlvl) + <(J + B-WsVa)] ] 
W2 {Da [J - Ws{Bl +pa~ PbH + B-WsVa] ~ J [B^ - Pb + W^rjl - Ta)} 

J [B-{Byvl + Blvl - Va) + Vl{pa - WWa + ^a)] 



Vb 



+ 



Da[J ~ Ws{Bl +pa- Pb)vl + B-WsVa] " J {Bl ^ Pb + W^rf^ - Ta) ' 



(4.25) 
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where it should be noted that equation ( I4.25t reduces to the corresponding hydrodynamical ex- 
pression in the Hmit of i? = [cf. equation (4.12) of Pons et al. (2000) or equation (3.13) of 
RezzoUa et al. (2003)]. Note also that using equations ( I4.16> and (I4.17> it is possible to ob- 
tain expressions for and in terms of the post-shock quantities p and J; the corresponding 
expressions are rather lengthy and uninspiring; for this reason we report them in AppendixlAl 

When all of the post-shock quantities are expressed as functions of only pt, and J (i.e. Vg), it is 
still necessary to express Vg as function of the post-shock pressure pt,. To do this we follow Pons 
et al. (2000) and use the original jump conditions ( 14. H . (^2) ™d (14. 3> to obtain 



where ((F)) = Fa + Fb and H = B^/J'^ - b"^ / is a shock invariant quantity (i.e. [[H]] = 0, 
I Anile 19S^ . Note that i3„ is not just the normal component of the magnetic field but, rather, the 
projection of h along n, i.e. 



S„ = 5^1^ = -lj+-^B- . (4.28) 



Equation i4.21l is also known as the Lichnerowicz adiabat, and represents the relativistic MHD 
counterpart of the Hugoniot adiabat. 

A couple of remarks should be made. Firstly, equations ( I4.26t - P.27t can be used for fast 
and slow shocks but not for an Alfven discontinuity. In this case, in fact, [[/i/p]] = 0, equa- 
tions (^26}-(|^^} are simple identities and the shock velocity Vg is trivially given by the local 
Alfven velocity . Secondly, for purely hydrodynamical shocks it is possible to find an ana- 
lytic expression for Vg as a function of the post-shock pressure [cf. equation (4.14) of Pons et al. 
(2000)]. In relativistic MHD, however, the corresponding analytic expression has not been found 
and equation (I4.26> needs to be solved numerically using a standard root-finding algorithm, but 
also increasing the computational costs considerably. To guarantee that we are using the right 
shock velocity, the root is searched in the approriate physical interval, i.e. \Vg\ G (|Va|, 1) for 
fast shocks and \Vg\ £ (|w^|, | Vyi|) for slow shocks. 



, (4.26) 



. P . 



. p 



= , (4.27) 



4.2. Solution across a rarefaction wave 



Rarefaction waves are self-similar solutions of the flow equations, i.e. equations in which all of 
the fluid quantities depend on x and t through the combination ^ = x/i. Using this as the inde- 
pendent variable, the set of partial differential MHD equations can be rewritten as the following 
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set of ordinary differential equations (ODEs) 



de 

d(r - b%°) d(S'^ - - Dt;^) 



de dC 

dC d^ 



de dC 

A{S^ - b%^) d{S^v=' - b^b^ 



df 



d^ 



dS^ 

d^ 

dB^ d(S^t;^ - B^w^) 



de 



(4.29) 
(4.30) 
(4.31) 
(4.32) 
(4.33) 
(4.34) 
(4.35) 
(4.36) 



Equation ( I4.29> can be further decomposed as 



c\ n n??"^ n7?^ ri7?^ 

{v^'-O^ + P [iv^ - OW'v^ + 1] _ + _ OpW'yy— + {v- - OpW^v^— = , 

(4.37) 

while combining equation ( I4.30> with equations (I4.31> -( l4.33> provides us with the relations 





-0— 


+ {l-^v 




-0— 

^' dC 


d? 




-0— 
dc 


^ de 



di ^ di ^ di di 



de 



de 



dC 



d^ 



dw" ^^^dp ^^^^d(6°6°) , ^^^^ d(b"6^) ^ ^d(b°b^) d(5^5^) 



dC 



dC 



dC 



dC 



Finally, rewriting the definition of the local sound speed 

hg dp / 

where s is the specific entropy, in terms of the self-similar variable 

dpg _ , odp 
d^ ~ '^'d^ ' 



0, 

(4.38) 

0, 

(4.39) 
. 

(4.40) 



(4.41) 



(4.42) 



and collecting the different terms in equations ( I4.29> -( |4.36> . we obtain the following system of 
seven ODEs in the seven variables ,v'^ , , B^, fully determining the solution across 
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a rarefaction wave 



ABy f B 



V 7] 



(4.43) 



dB^ 



(4.44) 



= (1 - v^O^ + + - 25^77 + W^v- - Oiw - V^)] ^ + 



2Bi 



= 



... , By[i-v^) 



dp 
d? 



du^ 



B'-{i-v'') 



dv^ _ (^"^ - 
"d^ 



de dC 



B^ 



B^By 



Bv 



(4.45) 



dfj " 



dBy , ..dS 



dC ' 
(4.46) 



= 



, dp 
d? 



,di;!' 



" B^ B^B^ 



dS?' 



Q = By'^-B^'±- + {v^-i)'±—., 
dC d^ d^ 

„,dw^ ^^dv"" , „ ^^di3^ 
d^ d^ d^ 



dS 



dC ' 
(4.47) 

(4.48) 



(4.49) 



The system of equations j4.43> - (l4.49t can be recast into a simple matrix form and non-trivial 
similarity solutions exist only if the determinant of the matrix of coefficients is zero. This condi- 
tion leads to a quartic equation in the self-similar variable ^ 



blcl-evlW^^{C-l)vlW'' + 



+ W^C (1 - W^) - {b^cl 



where 



b = 



and 



w 



C^ci + b\l-c 



e + 

(4.50) 
(4.51) 
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and whose roots coincide with the eigenvalues of the original system of equations \2Al - ^T^ . 
When — 0, equation (I4.50> reduces to a second-order equation whose roots provide the ve- 
locities of the left and right-going fast-waves. In the more general case when B'-^ ^ 0, however, 
the quartic cannot be recast as the product of two quadratic equations (as it is the case in New- 
tonian hydrodynamics) and the solution must be found numerically. The corresponding roots 
provide the velocities of the left and right-going slow and fast magnetosonic rarefaction waves, 
respectively. 

Using the appropriate root for ^, the system of ODEs ( l4.43> - fl.49> can be rewritten in terms 
of the total pressure to obtain a new system of six ODEs to be integrated from the value of 
pressure ahead the rarefaction to the one behind itf . The explicit expressions of these equations 
are rather lengthy and do not provide any important information; for this reason we report them 
in AppendixlEl 

4.3. Solution across an Alfven discontinuity 

The solution across Alfven discontinuities is found by imposing the continuity of p and p and 
then solving the system of equations (I4.15> -( l4.17> and ( I4.19> - (I4.20> . using Vg — Va, where Va = 
Vx+Bx/[W^{r|^^/w)] is the Alfven velocity for left (— ) and right (+) going waves, respectively. 
Since p and p are continuous across the Alfven discontinuity, a solution needs to be found only 
for the three components of v and for the tangential components of the magnetic field By and 
Bz- In general, and because no analytic solution was found, we solve the corresponding system 
of equations ( l4.15> -( l4TT7l . ( I4.19l l-( l4.20> numerically with a Newton-Raphson scheme. No major 
difficulties have been found in determining an accurate solution provided that the waves are all 
well separated and that a sufficiently accurate initial guess is provided (cf. solution in FigurelT^. 
For the latter we have used an approximate Riemann solver based on the Harten-Lax-van Leer- 
Einfeldt (HLLE) algorithm (Hart en et al. 19 83' Einfeldt 19881 and a moderate truncation error 
(i.e. using about 800 gridpoints for the tests reported here). However, considerable difficulties 
have been encountered if the waves are very close to each other This is the case, for instance, 
of test number 5 of Balsara (2001), in which the left-going Alfven discontinuity and the left- 
going slow rarefaction wave have very similar propagation velocities (cf. solution in Figure fTTl . 
The exact solution found in this case has a truncation error which is small, but larger that those 
reached in the other tests (cf. data in Table ll li . 

5. Tangential Magnetic Field Approach: "St-method" 

As done in Sect. |3 in what follows we describe in detail the approach referred to as the Bt- 
method, in which we calculate all of the variables in the Riemann fan using as unknowns the 
values of the tangential components of the magnetic field, i.e. B^ and B^. As mentioned in the 
Introduction, much of the inspiration in the development and use of this method comes from the 
corresponding approach developed by Ryu and Jones (1995) in nonrelativistic MHD. However, 
important differences are also present. 

In particular, in Newtonian MHD the problem can be solved using the norm of the tangential 

component of the magnetic field Bt = By + B^ and the rotation angle ifj = aTctan{B^ / B^) 
across Alfven discontinuities. This is because Bt is conserved across Alfven discontinuities and 
ip is constant across fast and slow-waves (see |Jeffrey 1966^ . As a result, the relevant system of 
equations is solved using as unknowns the values of Bt in regions R2-R3, R4-R5, R6-R7 of the 
Riemann fan in Figure|2]and the angle ip in regions R3-R6. At the contact discontinuity it is then 
necessary to solve a system of four equations, given by the continuity of v and of p, in the same 

f The number of equations to be solved reduces from seven to six because when using the total pressure 
as the self-similar variable one equation becomes then trivial, i.e. dp/ dp — 1. 
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four unknowns. This can be solved using root-finding techniques such as the Newton-Raphson 
method. Finally, when Bx — 0, the presence of only two fast waves and a tangential discontinuity 
makes the solution of the problem even simpler (see Ryu and Jones 1995 for details). 

In relativistic MHD, on the other hand, the value of Bt can be discontinuous across Alfven 
waves and the angle ip can vary across fast and slow-waves; it is then not possible to solve the 
system using the same method. Note also that the equations reported below both for shock and 
rarefactions waves are strictly valid only if B^ ^ and indeed should be used only in regions 
R4 and R5 of the Riemann fan shown in Figure|2] In these regions, only slow-waves are present 
and these do not appear when B^ = 0. 

5.1. Solution across a shock front 

To calculate the solution across a shock front within the Bt-method we start by considering the 
same system of equations in Section WA] but we solve equations (I4.1> -( l43t considering B^ 
and B^ as the unknown quantities. From equations ( I4.19> and ( I4.20> we express the post-shock 
values of and w^: 



V 

VI = 



1 

1 



(5.1) 
(5.2) 



Using now equation ( I4.13> to obtain the post-shock value of D 



and calculating the post-shock value of the total pressure using the invariance of hgBn, i-e. 
[[hgBn]] = (see Anile 1989), we can express all of the quantities as a function of the post- 
shock values of v^, B^, B^, and of the shock-velocity Vg. An analytic solution for the post- 
shock value of in terms of the other post-shock quantities was sought but not found, forcing 
to the numerical solution of one of the equations ( I4.15> - (I4.17> . Furthermore, in analogy with 
what done in the p-method, we calculate the value of the shock velocity by solving numerically 
equation (I4.26> . 

Finally, it may be useful to point out that the numerical solution of equation ( I4.26> is at times 
complicated by the existence of two acceptable roots in the interval of velocities in which the 
value of the slow shock velocity has to be found (i.e. between the value of and the value of the 
Alfven velocity). Because only one of these two roots will lead to a convergent exact solution, a 
careful selection needs to be made. The existence of these two roots could be related to a known 
problem in Newtonian MHD where the use of the tangential components of the magnetic field as 
the post-shock independent variables can lead to the presence of more than one solution (cf, for 
instance, Jeffrey and Taniuti 1964 1. This problem seems to be present also in relativistic MHD 
(Komissarov 2003 1, but it has not represented a serious drawback for the approach followed here. 
More work is needed to determine whether the use of the tangential components of the magnetic 
field as the post-shock independent variables is really optimal or whether different choices are 
preferable. 

5.2. Solution across a rarefaction wave 

To calculate the solution across a rarefaction wave within the Bf-method we use the same set 
of ODEs J4.29> - (I4.36> discussed in Section 14^21 with the only but important difference that we 
do not use ^ as self-similar variable but, rather, the norm of the tangential components of the 
magnetic field Bt. More specifically, we use equations J4.29t - (l4.31t together with equations 
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J4.35t -( l4.36t and substitute the derivative with respect to ^ with the one with respect to Bt- In 
addition to these equations, which provide a solution for variables p, p, w^, and if-, we express 
explicitly the relation between the norm and the tangential components in terms of the angle ip 

By = cos ipBt, (5.4) 

B^=sinV'Bt, (5.5) 
and rewrite them as ODEs having Bt as the self-similar variable 

cosV' , (5.6) 



dBy 



dBt 



sinV' . (5.7) 

Note that in deriving equations ( l5.6t -( l5T7t . an implicit assumption has been made: i.e. that the 
angle ip is constant across the rarefaction wave and thus that the tangential magnetic field does not 
rotate across the rarefaction wave. With the use of the supplementary equations (I5.6> -( I5^ . the 
resulting system of ODEs is complete and can be solved numerically using standard techniques 
for the solution of a system of coupled ODEs. In practice, the integration is started ahead of the 
rarefaction and is progressed toward the contact discontinuity, where Bt is given by the values of 
By and B^ chosen at the contact discontinuity. In all of the tests reported here (with the exception 
of test number 5 of Balsara 2001; see Section l6.2.1l for a discussion), the assumption — const, 
is valid. This is probably related to the choice of the initial conditions used in these tests and in 
particular to the fact that v\ = v\, B\ ~ B\, or v\ = B\ = 0, where A = {left, right), so that 
the initial states are essentially invariant after the exchange of y with z or the z components of v 
and B remain equal to zero in all the regions. 

It should be noted that also in relativistic hydrodynamics the velocity components tangential 
to a nonlinear wave can change their norm across the wave, in contrast with what happens in 
Newtonian hydrodynamics. Considering for simplicity the case for a shock wave in the limit of 
zero magnetic field, equations ( I4.16> -( l4.17> reduce to [[S'^/Z?]] = = [[S"^ /£)]], indicating that 
the ratio yy /v^ remains unchanged through shocks so that the tangential velocity 3-vector can 
change its norm but does not rotate. This property, which applies also across rarefaction waves, 
is not present across Newtonian nonlinear waves, in which the tangential 3-velocity vector does 
not rotate, nor changes its norm: [[w^]] = = [[w^]]. 

Although the condition ijj — const, is exact in nonrelativistic MHD, it may not be valid in 
relativistic regimes, where the tangential magnetic field is instead free to rotate across the slow 
rarefaction. In this case, a new strategy needs to be implemented and the simplest one consists 
of using the angle 'ip as the self-similar variable so that the system of equations ( I4.29> -( |4.36> can 
be expressed in terms of this new integration variable. In addition, the supplementary differential 
equation for one of the components of the tangential magnetic field can be obtained through the 
algebraic relation 

By = , (5.8) 

smtp 



and its derivative with respect to ^ 

dBy COS ip dB^ 



dip sin tp dtp sin ijj 



B"- . (5.9) 



The integration of the system of ODEs is done starting from the value of ip given by the ratio of 
the tangential components of the magnetic field ahead of the rarefaction wave, up to the value 
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given by the amplitudes of and at the contact discontinuity. Furthermore, as for the p- 
method, also within the i3t-method the values of the variable ^ are obtained from the quartic 
equation (I4.5Q> . 

A representative example of this effect is shown in Figure^] where we plot the exact solution 
of the generic Alfven test at time t = 1.5 (cf. Table0]for the initial conditions of this test). In 
particular, the left panel of Figure [O] shows the norm of the tangential magnetic field B^, while 
the right panel the angle ijj = arctan {B^ /B^). Note how both quantities vary across all the fast, 
slow and Alfven waves. 

6. Numerical Implementation and Representative Results 

Since the properties of the magnetic field components in the initial states lead to considerably 
different Riemann problems (cf. the two Riemann fans in Figures ^ and |3, we will discuss 
separately the numerical solution in the cases in which B^ = and B-^ ^ 0, emphasizing the 
properties of some of the most representative tests. 

6.1. Tangential Initial Magnetic Field: B'^ = 

As discussed in Section|3] when B^ = the Riemann problem consists of only two fast-waves 
and of a tangential discontinuity across which only and p are continuous {cf. Figure [Q. It 
should be noted that the condition of continuity of the total pressure across the tangential discon- 
tinuity does not necessarily extend also to the gas pressure and, indeed, the latter is in general 
discontinuous {cf. Figures |3l and |4}. In essence, the numerical solution of the Riemann problem 
when B'-^ = proceeds as follows: given the initial left and right states {i.e. regions Rl and 
R4 of Figure [0, we follow the procedure used in relativistic hydrodynamics and determine two 
unknown states as function of the common total pressure in regions R2 and R3 (p-method). The 
jump in the normal component of the velocity at the tangential discontinuity is then checked and 
a new guess for the total pressure found. This procedure is iterated until the solution is found 
with the desired accuracy. The numerical approach used is a combination of Newton-Raphson 
and bisection methods, starting from a value for the total pressure which is the average of the 
initial left and right states. Furthermore, to decide whether the wave considered is a shock or a 
rarefaction, we compare the values of the total pressure ahead of and behind the wave, solving the 
set of equations across a shock if the guessed value is larger than the total pressure ahead of the 
wave and thus in the initial state. We note that this procedure could be improved if an approach 
similar to the one discussed by Rezzolla et al. (2001, 2003) is implemented, which would exploit 
the values of the initial relative velocity to predict the wave-pattern produced. 

It is also worth noting that even though the numerical strategy discussed so far is very similar 
to the one used in relativistic hydrodynamics, the equations to be solved in MHD are much more 
complex and, more importantly, their computational cost markedly larger This is essentially 
because an analytic expression for the shock velocity was not found, so that the latter must be 
determined numerically. 

6.1.1. Representative Tests for B^ = 

Because initial states with a zero normal magnetic field lead to a Riemann problem that is 
comparatively much simpler to solve, an independent numerical code has been built for this case 
and it has been tested to reproduce known results in relativistic hydrodynamics, as well as a test 
proposed by Komissarov (1999) (this is referred to as the "shock-tube" test 2). We have also 
considered an additional, more generic shock-tube test in which all of the quantities in the initial 
states are nonzero and in which v ■ B ^ (this is referred to as the "generic shock-tube" test)f . 

f We note that a Riemann problem with — 0, but with v ■ B ^ cannot be solved with the exact 
solution recently proposed by Romero et al. (2005). 
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Test type 


P 


Pg 














Komissarov: Shock-'Ribe 2 (F = 4/3) 


















left state 


1.0 


30.0 


0.0 


0.0 


0.0 


0.0 


20.0 


0.0 


right state 


0.1 


1.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


Generic Shock-Tube (F = 5/3) 


















left state 


1.0 


0.01 


0.1 


0.3 


0.4 


0.0 


6.0 


2.0 


right state 


0.01 


5000 


0.5 


0.4 


0.3 


0.0 


5.0 


20.0 



Table 1. Initial conditions for the tests of the exact Riemann solver when the magnetic field has zero 

normal component, i.e. B^ — 0. 





P 


P 




yV 








Rl 


O.lOOOE+01 


0.2300E+03 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


0.2000E+02 


O.OOOOE+00 


R2 


0.2410E+00 


0.161 lE+02 


0.8497E+00 


O.OOOOE+00 


O.OOOOE+00 


0.9141E+01 


O.OOOOE+00 


R3 


0.6426E+00 


0.1611E+02 


0.8497E+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


R4 


O.lOOOE+00 


O.lOOOE+01 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 



Table 2. First significant digits for the exact solution of the test shock-tube 2 of Komissarov (1999) com- 
puted with an accuracy of 10^^^. The left column indicates the regions in which the solution is computed 
(cf Fig.0. 





P 


P 




yV 




B" 


B" 


Rl 

R2 
R3 
R4 


O.lOOOE+01 
0.1581E+01 
0.5489E-03 
O.lOOOE-01 


0.1819E+02 
0.4459E+02 
0.4459E+02 
0.5138E+04 


O.lOOOE+00 
-0.3073E+00 
-0.3073E+00 
0.5000E+00 


0.3000E+00 
0.3082E+00 
0.7488E+00 
0.4000E+00 


0.4000E+00 
0.2927E+00 
0.5556E+00 
0.3000E+00 


0.6000E+01 
0.9582E+01 
0.1023E+01 
0.5000E+01 


0.2000E+01 
0.3194E+01 
0.4092E+01 
0.2000E+02 



Table 3. The same as Table|2|but for the generic shock-tube test computed with an accuracy of 10 



Because the procedure for calculating the solution in this case is particularly simple and well 
tested from relativistic hydrodynamics, the algorithm employed has shown to be very robust and 
no failures were encountered in the calculation of any quantity. We list in Table^the set of initial 
conditions used in the tests solved, while we report in Tables|2]and|3lthe first significant digits for 
the exact solution of the same tests, reporting in all cases the accuracy obtained (which usually 
is < 10~^^). Finally, the full solutions in space of the various Riemann problems listed in Table 
[2 and for the quantities p, v^, pg, p, v^, v^, B^, and are shown in Figures |3] and |3 at the 
indicated representative times. 
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250 



200 



100 



50 



0.5 



-0.5 



0.5 



CQ 



-0.5 - 



Figure 3. Exact solution of the test shock-tube 2 of Komissarov (1999) at time t = 1.0. The solution 
composed of a left-going rarefaction wave, a tangential discontinuity and a right-going shock. 




Figure 4. Exact solution of the generic shock-tube test at time t = 1.0. The solution is composed of a 
left-going shock, a tangential discontinuity and a right-going rarefaction wave 
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6.2. Generic Initial Magnetic field: ^ 

As discussed in Section|3] when ^ the Riemann problem consists of seven different waves: 
two fast-waves, two slow-waves, two Alfven discontinuities and a central contact discontinuity 
across which only the density can be discontinuous (c/ Figure In essence, the numerical 
solution of the Riemann problem when B^ ^ proceeds as follows: starting from the initial left 
and right states {i.e. regions Rl and R8 of Figure|3, we compute the states after the fast-waves 
(regions R2 and R7), then we determine the jumps across the Alfven discontinuities (regions R3 
and R6) and finally we solve the equations for the slow-waves (regions R4 and R5). As a result 
of this sequence, the jump conditions in all the physical variables in the two states across the 
contact discontinuity are computed and if the solution obtained in this way does not reach the 
desired accuracy, the procedure is iterated. 

We also recall that when i?^ ^ 0, the numerical solution is found using a hybrid method 
which adopts different sets of equations according to the region in which the Riemann problem 
has to be solved. In particular, to compute the states after the fast-waves and across the Alfven 
discontinuities we use as unknown the total pressure (p-method; Section|4} and to discriminate 
between shocks and rarefaction waves we evaluate the jump in the total pressure in a way similar 
to the case when i?^ ~ 0. To compute the states after the slow-waves, on the other hand, we use 
the tangential components of the magnetic field (S^ -method; Section|5jl and to decide whether a 
wave is a shock or a rarefaction we evaluate the jump in the norm of the magnetic field bearing 
in mind that it must decrease across slow shocks and increase otherwise. Then at the contact 
discontinuity we compute the jumps in the total pressure and in the components of three-velocity 
and if they are above a certain accuracy we iterate by changing the values of the total pressure, 
used in regions R2-R3 and R6-R7, and of the tangential components of the magnetic field, used 
in regions R4-R5. 

It is worth underlining that the solution of the Riemann problem with generic initial states is 
considerably more demanding than when B^ = and not only because of the more numerous 
waves present. Indeed, the most severe difficulty is due to the fact that the set of equations to be 
solved becomes particularly stiff near the solution. A careful investigation of the several cases 
considered has in fact revealed that, in general, the functional behavior of the quantities whose 
roots are sought, changes very rapidly near the roots, stretching the ability of standards root- 
finding algorithms. As a result, it is not uncommon that the solution cannot be found if the 
iteration for the search of the root starts from a guess which is not sufficiently close to the exact 
solution. To avoid such failures and to provide a first guess which is reasonably accurate, we have 
used as a guide the solution provided by the HLLE approximate Riemann solverf . In practice, 
the approximate solution should be accurate to a few percent in the regions away from the waves, 
where the states are almost constant (very close to the waves the errors are of course much larger). 
Using this guess has proven to be sufficient to obtain a solution in all of the cases considered, 
but of course there is no guarantee that a solution will be straightforwardly found for all of the 
possible initial states. Our experience when the solution could not be immediately obtained is 
that an increase in the accuracy of the approximate Riemann solver is in general sufficient to 
yield a convergent and accurate solution. 

6.2.1. Representative Tests for B'^ 7^ 

Although the numerical code developed for the exact solution of the Riemann problem in 
relativistic MHD could in principle be used for generic initial data, we have used it in particular 
to calculate the exact solutions of those less trivial initial states that over the years have become 
standard references (e. g. IKomissarov 19991 or lBalsara 200 It . Table |3 collects the set of initial 

f Note that this is not necessary when = since in this case the solution can also be quite far from 
the exact one and yet the iterative scheme does not show problems in converging to it. 
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Test type 


P 
















Komissarov: Shock-Tube Test 1 (r = 4/3) 

left state 
right state 


1.0 
0.1 


1000.0 
1.0 


0.0 
0.0 


0.0 

0.0 


0.0 

0.0 


1.0 
1.0 


0.0 
0.0 


0.0 
0.0 


Komissarov: Collision Test (F = 4/3) 

left state 
right state 


1.0 
1.0 


1.0 
1.0 


5/\/26 
-5/^26 


0.0 
0.0 


0.0 
0.0 


10.0 
10.0 


10.0 
-10.0 


0.0 
0.0 


Balsara Test 1 (Brio & Wu) (F = 2) 

left state 
right state 


1.000 
0.125 


1.0 
0.1 


0.0 
0.0 


0.0 
0.0 


0.0 
0.0 


0.5 
0.5 


1.0 
-1.0 


0.0 
0.0 


Balsara Test 2 (F = 5/3) 
left state 
right state 


1.0 
1.0 


30.0 
1.0 


0.0 
0.0 


0.0 
0.0 


0.0 
0.0 


5.0 
5.0 


6.0 
0.7 


6.0 
0.7 


Balsara Test 3 (F = 5/3) 
left state 
right state 


1.0 
1.0 


1000.0 
0.1 


0.0 
0.0 


0.0 

0.0 


0.0 

0.0 


10.0 
10.0 


7.0 
0.7 


7.0 
0.7 


Balsara Test 4 (F = 5/3) 
left state 
right state 


1.0 
1.0 


0.1 
0.1 


0.999 
-0.999 


0.0 
0.0 


0.0 
0.0 


10.0 
10.0 


7.0 
-7.0 


7.0 
-7.0 


Balsara Test 5 (F = 5/3) 
left state 
right state 


1.08 
1.00 


0.95 
1.0 


0.40 
-0.45 


0.3 
-0.2 


0.2 
0.2 


2.0 
2.0 


0.3 
-0.7 


0.3 
0.5 


Generic Alfven Test (F = 5/3) 
left state 
right state 


1.0 
0.9 


5.0 
5.3 


0.0 
0.0 


0.3 
0.0 


0.4 
0.0 


1.0 
1.0 


6.0 
5.0 


2.0 
2.0 



Table 4. Initial conditions for the tests of the exact Riemann solver when the magnetic field has nonzero 

normal component, i.e. B^ 7^ 0. 



conditions used in the tests solved, while we report in Tablesl5t412lthe first significant digits for 
the exact solution of the same tests, reporting in all cases the accuracy obtained (which usually 
is ~ 10~^°). Finally, the full solutions in space of the various Riemann problems listed in Table 
|4]and for the quantities p, v^, pg, p, v^, v^, B^, and are shown in Figures I5lll2l at the 
indicated representative times. In addition. Figure [T3l offers a quantitative view of the changes 
in the tangential magnetic field Bf and of the rotation angle across the fast, slow and Alfven 
waves in the case of a generic Alfven test. 

In all of the tests reported in Table |4] the HLLE solver with about 800 gridpoints was able 
to track rather well the exact solution in all of its waves. The only exception to this has been 
test number 1 of Balsara (2001) which represents the relativistic version of the test proposed 
by Brio-Wu (1988) in Newtonian hydrodynamics ( van Putten 1993 1. The approximate numerical 
solution of this test, in fact, shows the development of a left-going slow compound-wave, that is 
a wave composed by a slow shock adjacent to a slow rarefaction. Since we assume that a slow 
or fast-wave can either be a pure rarefaction or a pure shock, compound structures of this type 
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P 


P 




V 


z 

V 






Rl 


O.lOOOE+01 


O.lOOlE+04 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


R2 


0.6984E-01 


0.2927E+02 


0.9115E+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


R3 


0.6984E-01 


0.2927E+02 


0.9115E+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


R4 


0.6984E-01 


0.2927E+02 


0.9115E+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


R5 


0.8846E+00 


0.2927E+02 


0.9115E+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


R6 


0.8846E+00 


0.2927E+02 


0.9115E+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


R7 


0.8846E+00 


0.2927E+02 


0.9115E+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


R8 


O.lOOOE+00 


0.1500E+01 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 



Table 5. First significant digits for the exact solution of the test shock-tube 1 of Komissarov (1999) com- 
puted with an accuracy of 10~^". The left column indicates the regions in which the solution is computed 
{cf. Fig.H}. 



cannot be found by construction and thus are not present in the exact solution found (cf. Table0 
and Figure Q. We remark that it is not yet clear whether compound waves have to be consid- 
ered acceptable physical solutions of the ideal MHD equations and a debate on this is still on- 
going (see, for instance, Myong & Roe IQQTa'/Myong & Roe 1997b Falle & Komissarov 206T\ 
ITorrilhon 2003allTorrilhon 2003b. Torrilhon & Balsara 2004J . We here prefer to adopt the same 
standpoint of Ryu and Jones (1995) in the development of their exact Riemann solver in non- 
relativistic magnetohydrodynamics and not comment further on this until a commonly accepted 
view has emerged. 

Another test which deserves a special comment is test number 5 of Balsara (2001 ), in which the 
left-going Alfven discontinuity and the left-going slow rarefaction wave have very similar prop- 
agation velocities. Indeed they are so close to each other that not even the HLLE approximate 
Riemann solver with 40000 gridpoints was able to capture the precise location of the discontinu- 
ity. As a consequence, the initial guess for the jumps across the left-going Alfven discontinuity 
was sufficiently good to yield a convergent solution, but not good enough to provide an exact 
solution with a truncation error comparable with the one reached in all of the other tests (cf. data 
in Table ll II . In addition, another distinctive feature of this test and which has not been found in 
any of the others, is the rotation of the angle i/; across the left-going slow rarefaction. To handle 
this we have followed the procedure discussed in Section ls!2l and used equation i5.9\ to compute 
the changes in the tangential magnetic field across the rarefaction wave. 

7. Conclusions 

We have presented the procedure for the solution of the exact Riemann problem in special rel- 
ativistic MHD. Special care has been paid in treating both degenerate initial states (i.e. with zero 
normal magnetic field) leading to a set of only three waves analogous to the ones in relativistic 
hydrodynamics, as well as generic initial states (i.e. with nonzero normal magnetic field) leading 
to the full set of seven MHD waves. 

The approach discussed for the numerical solution of the exact Riemann problem reflects this 
distinction and different sets of equations are used according to the values of the normal magnetic 
field. In particular, when — 0, all of the equations needed for the solution of the Riemann 
problem are written as a function of the total pressure, thus following a procedure which is 
logically equivalent to the one adopted in relativistic hydrodynamics (we have referred to this as 
to the p-method). When ^ 0, on the other hand, an hybrid approach is adopted in which the 
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Figure 5. Exact solution of the test shock-tube 1 of Komissarov (1999) at time t = 1.0. The solution is 
composed of a left-going fast rarefaction, of a contact discontinuity and of a right-going fast shock. 
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Figure 6. Exact solution of the collision test of Komissarov (1999) at time t = 1.22. The solution 
is composed of a left-going fast shock, of a left-going slow shock, a right-going slow shock and of a 
right-going fast shock. 
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Figure 7. Exact solution of the test number 1 of Balsara (2001) at time t = 0.4 and which represents 
the relativistic version of the Brio-Wu test (1988). The solution is composed of a left-going fast rarefaction, 
of a left-going slow shock, of a contact discontinuity, of a right-going slow shock and of a right-going fast 
rarefaction. Note the absence of a slow compound-wave which cannot be found by construction in our exact 
solver, but that appears in the solution of the HLLE approximate Riemann solver (not shown). 




Figure 8. Exact solution of the test number 2 of Balsara (2001) at time t = 0.4. The solution is composed 
of two left-going fast and slow rarefactions, of a contact discontinuity and of two right-going fast and slow 
shocks. 
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Figure 9. Exact solution of the test number 3 of Balsara (2001) at time t = 0.4. The solution is composed 
of two left-going fast and slow rarefactions, of a contact discontinuity and of two right-going fast and slow 
shocks. 
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Figure 10. Exact solution of the test number 4 of Balsara (2001) at time t = 0.4. The solution is 
composed of two left-going fast and slow shocks and of two right-going fast and slow shocks. 
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Figure 11. Exact solution of the test number 5 of Balsara (2001) at time t = 0.55. The solution is 
composed of a left-going fast shock, of a left-going Alfven discontinuity, of a left-going slow rarefaction, 
of a contact discontinuity, of a right-going slow shock, of a right-going Alfven discontinuity and of a 
right-going fast shock. Note that the accuracy in this test is only rather low: 3 x 10~*. 
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Figure 12. Exact solution of the generic Alfven test at time t — 1.5. The solution is composed of a 
left-going fast rarefaction, of a left-going Alfven discontinuity, of a left-going slow shock, of a contact 
discontinuity, of a right-going slow shock, of a right-going Alfven discontinuity and of a right-going fast 
shock. 
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Figure 13. Exact solution of the generic Alfven test at time t — 1.5. The left panel shows the norm of 
the tangential magnetic field Bt, while the right panel the angle tp = arctan (B^/B"). Note that both 
quantities vary across all the fast, slow and Alfven waves as a result of a relativistic effect. 





P 


P 




yV 








Rl 


O.lOOOE+01 


0.5292E+02 


0.9806E+00 


O.OOOOE+00 


O.OOOOE+00 


O.lOOOE+02 


O.OOOOE+00 


R2 


0.6331E+01 


0.257 lE+03 


0.4380E+00 


0.4069E+00 


O.OOOOE+00 


0.1960E+02 


O.OOOOE+00 


R3 


0.6331E+01 


0.257 lE+03 


0.4380E+00 


0.4069E+00 


O.OOOOE+00 


0.1960E+02 


O.OOOOE+OO 


R4 


0.2742E+02 


0.2819E+03 


0.2453E-07 


-0.6811E+00 


O.OOOOE+00 


0.2250E-06 


O.OOOOE+OO 


R5 


0.2742E+02 


0.2819E+03 


-0.2810E-07 


-0.6811E+00 


O.OOOOE+00 


0.2250E-06 


O.OOOOE+OO 


R6 


0.6331E+01 


0.257 lE+03 


-0.4380E+00 


0.4069E+00 


O.OOOOE+00 


-0.1960E+02 


O.OOOOE+00 


R7 


0.6331E+01 


0.257 lE+03 


-0.4380E+00 


0.4069E+00 


O.OOOOE+00 


-0.1960E+02 


O.OOOOE+00 


R8 


O.lOOOE+01 


0.5292E+02 


-0.9806E+00 


O.OOOOE+00 


O.OOOOE+00 


-O.lOOOE+02 


O.OOOOE+00 



Table 6. The same as Table|5]but for the exact solution of the test Collision of Komissarov (1999) 

computed with an accuracy of 10^''. 





P 


P 


X 

V 


V 

V 


z 

V 




B" 


Rl 


O.lOOOE+01 


0.1625E+01 


O.OOOOE+OO 


O.OOOOE+OO 


O.OOOOE+OO 


O.IOOOE+Ol 


O.OOOOE+00 


R2 


0.6257E+00 


0.6989E+00 


0.3742E+00 


-0.3561E-01 


O.OOOOE+OO 


0.6594E+00 


O.OOOOE+00 


R3 


0.6257E+00 


0.6989E+00 


0.3742E+00 


-0.3561E-01 


O.OOOOE+OO 


0.6594E+00 


O.OOOOE+OO 


R4 


0.7092E+00 


0.7062E+00 


0.2555E+00 


-0.6804E+00 


O.OOOOE+00 


-0.4285E+00 


O.OOOOE+OO 


R5 


0.2695E+00 


0.7062E+00 


0.2555E+00 


-0.6804E+00 


O.OOOOE+00 


-0.4285E+00 


O.OOOOE+00 


R6 


0.1223E+00 


0.6976E+00 


-0.2080E-01 


-0.3460E-02 


O.OOOOE+OO 


-0.9769E+00 


O.OOOOE+00 


R7 


0.1223E+00 


0.6976E+00 


-0.2080E-01 


-0.3460E-02 


O.OOOOE+OO 


-0.9769E+00 


O.OOOOE+OO 


R8 


0.1250E+00 


0.7250E+00 


O.OOOOE+OO 


O.OOOOE+OO 


O.OOOOE+OO 


-O.lOOOE+01 


O.OOOOE+OO 



Table 7. The same as Table|5|but for the exact solution of the test number 1 of Balsara (2001) computed 
with an accuracy of 10"^". This test represents the relativistic version of the test proposed by Brio & Wu 
(1988) 
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P 


P 




V 


z 

V 






Rl 


O.lOOOE+01 


0.7850E+02 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


0.6000E+01 


0.6000E+01 


R2 


0.4300E+00 


0.2321E+02 


0.6344E+00 


-0.998 lE-01 


-0.998 lE-01 


0.3045E+01 


0.3045E+01 


R3 


0.4300E+00 


0.2321E+02 


0.6344E+00 


-0.998 lE-01 


-0.998 lE-01 


0.3045E+01 


0.3045E+01 


R4 


0.3830E+00 


0.2284E+02 


0.6770E+00 


-0.5566E-01 


-0.5566E-01 


0.3205E+01 


0.3205E+01 


R5 


0.2828E+01 


0.2284E+02 


0.6770E+00 


-0.5566E-01 


-0.5566E-01 


0.3205E+01 


0.3205E+01 


R6 


0.1582E+01 


0.2072E+02 


0.4688E+00 


-0.2538E+00 


-0.2538E+00 


0.3971E+01 


0.397 lE+01 


R7 


0.1582E+01 


0.2072E+02 


0.4688E+00 


-0.2538E+00 


-0.2538E+00 


0.3971E+01 


0.397 lE+01 


R8 


O.lOOOE+01 


0.1399E+02 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


0.7000E+00 


0.7000E+00 



Table 8. The same as Table|5|but for the exact solution of the test number 2 of Balsara (2001) computed 

with an accuracy of 10^^". 





P 


P 








B^ 




Rl 


O.lOOOE+01 


0.1099E+04 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


0.7000E+01 


0.7000E+01 


R2 


0.1381E+00 


0.8604E+02 


0.9246E+00 


-0.3513E-01 


-0.3513E-01 


0.2238E+01 


0.2238E+01 


R3 


0.1381E+00 


0.8604E+02 


0.9246E+00 


-0.3513E-01 


-0.3513E-01 


0.2238E+01 


0.2238E+01 


R4 


0.9798E-01 


0.7653E+02 


0.9529E+00 


0.4366E-01 


0.4366E-01 


0.4670E+01 


0.4670E+01 


R5 


O.lOlOE+02 


0.7653E+02 


0.9529E+00 


0.4366E-01 


0.4366E-01 


0.4670E+01 


0.4670E+01 


R6 


0.1218E+01 


0.6363E+02 


0.4670E+00 


-0.4270E+00 


-0.4270E+00 


0.9408E+01 


0.9408E+01 


R7 


0.1218E+01 


0.6363E+02 


0.4670E+00 


-0.4270E+00 


-0.4270E+00 


0.9408E+01 


0.9408E+01 


R8 


O.lOOOE+01 


0.5059E+02 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


0.7000E+00 


0.7000E+00 



Table 9. The same as Table|5|but for the exact solution of the test number 3 of Balsara (2001) computed 

with an accuracy of 10^^". 



solution across fast-waves and Alfven discontinuities is still computed using the p-method, but 
the one across slow-waves and the contact discontinuity is computed using equations which are 
written in terms of the tangential components of the magnetic field (we have referred to this as 
to the i?t-method). The use of a combined approach for the general case of ^ has turned 
out to be crucial for a successful solution of the problem. 

Because of its generality, the solution presented here could serve as a useful if not indispens- 
able test for those numerical codes that solve the MHD equations in relativistic regimes. As 
the astronomical observations become increasingly more accurate, such numerical codes will 
become increasingly more important to explain and describe in detail the complex physics of 
astrophysical compact objects. 

As a final remark we note that despite the considerable improvements in the performance 
of modern computers, the exact solution of the Riemann problem at each grid interface is still 
computationally too expensive to be used routinely in sophisticated multidimensional numerical 
codes solving the equations of relativistic hydrodynamics or MHD in either stationary or dynam- 
ical spacetimes (see, for instance, Baiotti et al. 2005, Duez et al. 2005). While a numerical code 
based on exact Riemann solvers may represent at least in principle the most accurate approach 
to the solution of the hydrodynamics and MHD equations, considerable work is still required to 
make it competitive with less accurate but more computationally efficient methods. A first step 
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P 


P 


X 

V 


V 

V 


z 

V 


B^ 


B" 


Rl 


O.lOOOE+01 


0.5020E+02 


0.9990E+00 
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0.7000E+01 


0.7000E+01 


R2 


0.5175E+02 


0.1184E+04 


0.4408E-01 


0.3263E-01 


0.3263E-01 


0.1668E+02 


0.1668E+02 


R3 


0.5175E+02 


0.1184E+04 


0.4408E-01 


0.3263E-01 


0.3263E-01 


0.1668E+02 


0.1668E+02 


R4 


0.6148E+02 


0.1188E+04 


0.1086E-07 


-0.2877E+00 


-0.2877E+00 


0.8042E-09 


0.8036E-09 


R5 


0.6148E+02 


0.1188E+04 


-0.1089E-07 


-0.2877E+00 


-0.2877E+00 


0.8042E-09 


0.8036E-09 


R6 


0.5175E+02 


0.1184E+04 


-0.4408E-01 


0.3263E-01 


0.3263E-01 


-0.1668E+02 


-0.1668E+02 


R7 


0.5175E+02 


0.1184E+04 


-0.4408E-01 


0.3263E-01 


0.3263E-01 


-0.1668E+02 


-0.1668E+02 


R8 


O.lOOOE+01 


0.5020E+02 


-0.9990E+00 


O.OOOOE+00 


O.OOOOE+00 


-0.7000E+01 


-0.7000E+01 



Table 10. The same as Table|5|but for the exact solution of the test number 4 of Balsara (2001) computed 

with an accuracy of 10^^. 





P 


P 












Rl 


0.1080E+01 


0.2885E+01 


0.4000E+00 


0.3000E+00 


0.2000E+00 


0.3000E+00 


0.3000E+00 


R2 


0.2447E+01 


0.5908E+01 


-0.1331E+00 


0.2111E+00 


0.1751E+00 


0.2662E+00 


0.5076E+00 


R3 


0.2447E+01 


0.5908E+01 


-0.1215E+00 


0.1264E+00 


0.1158E+00 


-0.1182E+00 


0.2302E+00 


R4 


0.2050E+01 


0.5616E+01 


-0.4547E-01 


-0.1463E+00 


0.2146E+00 


-0.1175E+01 


0.5852E+00 


R5 


0.1884E+01 


0.5616E+01 


-0.4543E-01 


-0.1462E+00 


0.2149E+00 


-0.1175E+01 


0.5850E+00 


R6 


0.1642E+01 


0.5488E+01 


-0.1129E+00 


-0.4606E-01 


0.1601E+00 


-0.1429E+01 


0.7320E+00 


R7 


0.1642E+01 


0.5488E+01 


-0.1155E+00 


-0.8536E-01 


0.1027E+00 


-0.1272E+01 


0.9468E+00 


R8 


O.lOOOE+01 


0.2918E+01 


-0.4500E+00 


-0.2000E+00 


0.2000E+00 


-0.7000E+00 


0.5000E+00 



Table 1 1 . The same as Table|5|but for the exact solution of the test number 5 of Balsara (2001) computed 

with an accuracy of 3 x 10"*. 



in this direction would be, for instance, the search for an analytic solution for the shock velocity 
and this will be the subject of future work. Another important problem deserving equal attention 
is that of the uniqueness of the solution. While a global consensus on this issue still needs to be 
reached, it will remain essential in order to construct a complete and consistent picture of the 
exact solution of the Riemann problem in relativistic MHD. 

The numerical codes computing the exact solution both when — and when ^ are 
available from the authors upon request. Users of the codes can give credit by mentioning the 
source and citing this paper. 

It is a pleasure to thank Jose M-. Marti, Jose A. Pons and OUndo Zanotti for useful discussions 
and comments. 
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Rl 


O.lOOOE+01 


0.2376E+02 


O.OOOOE+00 


0.3000E+00 


0.4000E+00 


0.6000E+01 


0.2000E+01 


R2 


0.9219E+00 


0.2083E+02 


0.6232E-01 


0.3050E+00 


0.4193E+00 


0.5622E+01 


0.1892E+01 


R3 


0.9219E+00 


0.2083E+02 


0.7109E-01 


0.3669E+00 


0.2429E+00 


0.5691E+01 


0.8502E+00 


R4 


0.1263E+01 


0.2087E+02 


0.3886E-01 


0.1147E+00 


0.2054E+00 


0.5130E+01 


0.7680E+00 


R5 


0.1099E+01 


0.2087E+02 


0.3886E-01 


0.1147E+00 


0.2054E+00 


0.5130E+01 


0.7680E+00 


R6 


0.9130E+00 


0.2085E+02 


0.1607E-01 


-0.5009E-01 


0.1813E+00 


0.5505E+01 


0.8195E+00 


R7 


0.9130E+00 


0.2085E+02 


0.1341E-01 


-0.6599E-03 


-0.2640E-03 


0.5073E+01 


0.2029E+01 


R8 


0.9000E+00 


0.2030E+02 


O.OOOOE+00 


O.OOOOE+00 


O.OOOOE+00 


0.5000E+01 


0.2000E+01 



Table 12. The same as Table|5|but for the exact solution of the generic Alfven test computed with an 

accuracy of 10"^". 



Appendix A. 

We here report the expressions for the tangential components of the velocity behind the shock 
(i.e. vf, v^) when expressed as function of post-shockp and J. First, we consider as function 
of Pb, J and v^. 



vl{Da+Pa +Ta- vX)iDa + Pb + Ta - 
By(B.,{pa + Ta)vl + Blipb + Ta)vl ~ r],{pb + T,) - 

Dairia - B^vl - Blvl) + rilWl ~ ^l{B,vl + B^vDWl) - 

J^DaWl[{2Blj^aVl + Blvl[{Da+Pa + ~ 2Bll^a] + 

<iPa - PbMDa +Pa+ ra)w^ + By{B'^v', - 2r/,)] - B,{7^a{3Da + 2pa + Pb + 3t,)«^ - 
B^^Pa-pbHv^b + BliiPa + ra){vf - 1) - 377^ + 2B^,7JaV^, + [pb + T<,)<«,^ - MDa]})W^, 

BrBliDa +Pb + Ta) + vl [2(P6 - PaKvy - B^v^vl + B,{ir^avl - MBl)]W^^ 



Blvf,iBy + 2r,avyw!)-Btvy + 

B.vlipa - Pb)[2By + {2^,vy + Blvlvl + BlvlvDW^,] + 

Bl(By{Blvl - 3fja + 2MvaW!) + vy[Da+Pa + r„ - + {pb - PaK^DW! 
B,Dl(BlvlvlWl - B^{pa-pb)vlvlvlvlWt - Q{Bl + r^avlWl) - 
BllBl + [t^^vI - MBy)W^^])w^] , 
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where we have defined M = (1 — v^^ — v^v^), Q = W^{pa — Pb){v%^ — 1) and 



iVa{Da +Pb + Ta) + 'irjiWl + Blvl(Da - Pa + 2pb + Ta ^ vlW^)] + 
2DaJ{Pa - Pb)vl{2na - Blvl)Ws + D^r?a(Pa " Pb){vf - l)wf) - 

BtDlW^iJ + DaV^W,) - WMDa+Pb + Ta- rfjVl) + 

DaiPb - Pa)vlWs][j\D, +Pb + Ta- + 2DaJ{pb - Pa)vlWs - 

Bl^J^ + Dlipt, - pa){vf - 1)W^] + BlDM2naJ^Wl + 
2D^rjaJvlWlWs + Dlirja - -B^^)^^,^] + B^ [j^W^{Da + P6 + t„ - vX) + 

DaJ^vlWliDa - 2p, + 3P6 + Ta- vX)W, + 

DIj(d^ +Ph- Bl^ + r„ + [2Bl'navl - - (Pa - Pb){Svf - 1 + v^^W^) 

Dlv2vfW^{p,-Pa)W^]} . 



We next consider the expression of as function of post-shock p and J 



= 'Wj-^"^' +Pa + Ta - vX){Da ' B^' - By' + 

Pb + Ta- rilWl) + Bl + Ta){B^vl + Blvl) - rjaiPb + Ta) + 

DaiB-vl + Blvl rja) + riw! vl{B^%: + Blvl)W^)] + 
DaJ'w! \B^B^,{Da+Pb + To) + (B^Bl[2Blriavl - 3vl + 



Daivf + Vf - 1) + {pa + Ta){vf + vf 1)] - 

Bl{Da-rPa + Ta)vl\vl - B-\:[{Da+Pa + " ^KVa] ' 

2{pa-pbK[{Da+Pa+TaK " B'^W^ + 



1) - {B^\l - SB^f^a + 2pavl - 2pbvl + B^Blvl)vl]Wt 



7^l(B-Bl{vf 

DlJWl (2B-{pa - Pb)vl{Bl + navlWl) - B-\l + B^\l{Bl + 2r,avlWl) 

B^\vl[Da - Bl' +Pa+Ta- Vai^Va - ^Blv^W^] + 

Bl[Blvl - 3r7„ - 2^a{vf + vl' - l)Wl\} + 

{Pa-Pb){vf - l)Wi[{Da+Pa+Ta)vl - % (i?^ + JJa^^^Ol) + 

B^Dl [bIvIvIWI - {pa - Pb){vf - 1)W^{B: + VaV^W^) - 
Bl + mvf + vf - 1) + - Blvl)vmt)\ Wl\ , 



+ 
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where 

BlvliD, -pa + 2p, + Ta- vlW^) + Blvl{Da - Pa + 2pb + r„ - ^IWl)] + 

2DaJ{pa - Pb)v:{27Ja " i?^^ " B^OW^. + Ol^JaiPa " Pb)ivf " 1)W^} - 
BtDlW^J + DaVlWs) - Wl[J(Da +Pb+Ta- " 
Daipa - Pb)vlWs][j\Da - Bl^ - Bf + Pb + Ta - vlW^) - 
2DaJ{Pa - PbKWs - DliPa - Pb){vf - 1)W^] + 

BlDaWs[2riaJ^Wl + 2DaVaJv:W^W, + Dl{r]a - Blvl - Blvl)Wl\ + 

Bl (j^WKPa + P6 + T„ - y)lWl) + BaJ\lWl(Pa - 2pa + Spb + T„ - vlW^)W, + 

DlJ{Da - Bf - Bf +Pb + Ta + [2Va{Blvl + Blv^) - 3ril - 

{Pa - Pb){3vf + Vf + Vf - 1)]W^}W!) } . 



Appendix B. 

The explicit form for the system of ODEs to be solved numerically to determine the solution 
across a rarefaction wave within the p-method is given by the following set of equations in which 
the total pressure p plays the role of the self-similar variable 



dp 



dvz 
dp 



= -plW\, + \ ^ - pW\-^ - pW^v,-^ , (Bl) 

dp V -t,J dp dp dp 



^ = R{{ph,W' + Bl){^ - v.){v.^ - 1) + Bl^^^j-^^ + BUivl + vl) + 

BMi.'' -^)-B.Vx{l-2v^^ + e)]} , (B2) 



R{ 2B,Vy{Tj - B,v,)^ - Blvyiii + v^) + 



Vy[Bl + W^jf - w)]{v, - 0^ + Blvy{u,X - 1) + ByB,v,{e - 1) + 

{vl + vl - 1) + {v, - 2v^vi)^ +{i+vi- vDe - v.e ^ 

{v. - 



B^By '-^-^ — '—^ — ^ " — '-^ — — \, (B 3) 



= i?|2B,(77 - ByVy)v,X - Blv,i{v, + e) + 

v,[Bl + 14^2(^2 _ _ ^ ByB,Vy{e - 1) + v,Bl{v^£, - 1) + 

R „ {vl + - 1) + K - 2v,vi)£, ^{i-vi+ vDe - v.e \ 

BxBz -T > , (.B4) 

{Vx - 



dBy _ W^By ~ ByV^C + B^Vy^) 



dp Bl + 2B,7^W^{v^, - + 1^4(^2 _ ^)(^^ _ ^)2 

d^. W^{B,- B,v^^ + B,v,i) 

dp Bl + 2B,r]W^v, - + - w){vx - 



(B5) 
(B6) 
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where we have defined 



with 



yf=«, + — ^— (B8) 



being the Alfven velocities in the two directions. Note that the set of ODEs has a singular point if 
the characteristic velocity of the slow or fast magnetosonic waves is equal to the Alfven velocity 
[cf. eq. ( IB 7H and cannot be solved in this case without a proper regularization. This procedure 
is not included in the numerical code made available upon request. 
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